NASA Technical Memorandum 86382 

. j 

■ j 

] NASA-TM-86382 19850025225 

i 

i 

- I 

“ ' , ■ A .. 

Selecting Step Sizes 

in Sensitivity Analysis 

by Finite Differences UlftW/v O^'L 


Jocelyn Iott, Raphael T. Haftka, 
and Howard M. Adelman 

AUGUST 1985 


NASA 





NASA Technical Memorandum 86382 


Selecting Step Sizes 
in Sensitivity Analysis 
by Finite Differences 


Jocelyn Iott and Raphael T. Haftka 

Virginia Polytechnic Institute and State University 

Blacksburg, Virginia 

Howard M. Adelman 
Langley Research Center 
Hampton, Virginia 


NASA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Branch 


1985 



Summary 

This paper deals with methods for obtaining near- 
optimum step sizes for finite difference approxima- 
tions to first derivatives with particular application 
to sensitivity analysis. A technique denoted the fi- 
nite difference (FD) algorithm, previously described 
in the literature, is reviewed and applied in a tuto- 
rial manner to the derivative of a sine function. The 
original FD algorithm is applicable to one derivative 
at a time. This paper describes an extension of the 
FD method to the calculation of several derivatives 
of matrix equations. 

Both the original and extended FD algorithms 
were applied to sensitivity analysis for a data-fitting 
problem in which a polynomial is passed through 
data points, and the derivatives of the coefficients of 
the polynomial with respect to uncertainities in the 
data were calculated. The methods also were applied 
to sensitivity analysis of the structural response of 
a finite-element-modeled swept wing. In a previous 
paper, this sensitivity analysis of the swept wing 
required a time-consuming trial-and-error effort to 
obtain a suitable step size, but it proved to be a 
routine application for the FD algorithm herein. 

Introduction 

One of the key areas in the development of meth- 
ods for optimization of engineering systems is sensi- 
tivity analysis — the calculation of derivatives of sys- 
tem response with respect to design parameters. In 
gradient-based optimization methods, derivatives are 
used to determine design changes to move toward 
an optimum design. Sensitivity derivatives also are 
used for rapid assessment of the effect of a change in 
a design parameter on the response of an engineer- 
ing system. In these typical applications, especially 
when a large number of repetitive calculations are re- 
quired, the finite difference method is often a useful 
technique for computing the derivatives. 

One problem which arises in connection with the 
finite difference method is the choice of step size. 
The step size can contribute to two types of errors 
in the finite difference approximation — truncation er- 
ror and condition error. Truncation error is the dif- 
ference between the exact value and the computed 
value of a perturbed function if it is assumed that 
there is no loss of numerical precision in the calcu- 
lation. Truncation error is generally represented by 
neglected terms in the Taylor series representation 
of a perturbed function. Truncation errors are in- 
creased by overly large step sizes. Condition errors 
are associated with numerical noise and are caused 
by loss of numerical precision. These errors may re- 
sult from computer round-off error or the operation 


of subtracting large numbers which are very nearly 
equal. Condition errors in finite difference derivatives 
generally increase with decreasing step sizes. 

In the absence of methods for selecting a suitable 
step size, a trial-and-error approach, which requires 
many function evaluations for each derivative, must 
be used (ref. 1). Automatically calculating an op- 
timum step size can decrease the time and cost of 
the calculations. One such method is described in 
reference 1. This method, referred to as the finite 
difference (FD) algorithm, is formulated on the basis 
that the finite difference error is minimized when the 
truncation error is equal to the condition error. The 
FD algorithm has previously been demonstrated for 
single derivatives and proved to be effective in pro- 
ducing step sizes that gave good approximations to 
the exact derivatives (refs. 1 and 2). The need to 
apply the method to functions governed by matrix 
equations has led to an extension of the FD algo- 
rithm and is the main subject of this paper. 

In this paper, the FD algorithm is applied to 
three examples: (1) the sine function, (2) a matrix 
equation (ref. 3, pp. 12-13) involving data fitting by 
a polynomial in which derivatives of coefficients with 
respect to uncertainties in the data are calculated, 
and (3) a finite element structural analysis of a 
swept wing. The polynomial and wing problems 
are used to test both the original and the extended 
versions of the FD algorithm. The finite element 
model of the swept wing is included in this study 
because in a previous sensitivity study, described in 
reference 4, the uncertainty of step size became a 
major concern, and a succession of trials was required 
to determine an acceptable step size. This paper 
contains a summary of the original FD algorithm, 
the development of the extended FD algorithm, and 
the results of applying the algorithms to the above 
problems. 

Symbols 

a,- coefficient of ith term of iVth-order 

polynomial 

C(h) condition error in approximating first 

derivative by finite difference method 

c(<f>) condition error in approximating 

second derivative by finite difference 
method 

e bound on total error in approximating 

derivatives by finite difference method 

F,F force vectors 

/ exact value of a function 



/' exact first derivative of a function 

f" exact second derivative of a function 

fU bound on second derivative 

h step size 

h op t optimum step size 

h s step size used for estimating second 

derivative 

K stiffness matrix 

N degree of polynomial 

T(h) truncation error in approximating first 
derivative by finite difference method 

U, U exact and computed displacement 
vectors, respectively 

x independent variable 

X{, t/i ith data point in Vandermonde 

problem 

e A bound on absolute error in computed 

function values 

$ finite difference approximation to 

second derivative 

Subscripts: 

max maximum 

min minimum 

Summary of FD Algorithm 
Truncation and Condition Errors 

The FD algorithm (ref. 2) is used to calculate an 
optimum step size for a first derivative approximated 
by forward differences 

/'M = /(x+ ' , ^ /(x) (i) 

Whenever finite difference formulas are used to ap- 
proximate derivatives, there are two sources of error: 
truncation error and condition error. The trunca- 
tion error T(h) is a result of the neglected terms in 
the Taylor series expansion of the perturbed func- 
tion and is represented in references 1 and 2 as being 
proportional to the step size as follows: 

T{h) = ^f"(x + (h) 0<<r<l (2) 

The condition error is the difference between the nu- 
merical evaluation of a function and the exact value. 


One contribution to the condition error is round-off 
error in evaluating equation (1), which is compara- 
tively small for most mainframe computers and is, 
therefore, considered negligible. However, if f(x) is 
computed by a lengthy or ill-conditioned numerical 
process, the round-off contribution to the condition 
error can be substantial. Additionally, condition er- 
rors may result if f(x) is calculated by an iterative 
process which is terminated early. The condition er- 
ror C{h ) can be conservatively estimated as being 
inversely proportional to the step size. The following 
expression, used herein, is given in reference 2: 

C(h) = \e A (3) 

Determination of e A 

In equation (3), e A is a bound on the absolute 
error in the computed function /. Obtaining an ap- 
propriate value of e A for use in equation (3) gen- 
erally requires some knowledge of the problem be- 
ing solved. A straightforward (but not always con- 
venient) method of obtaining e A is to perform the 
calculation of / in single and double precision and 
subtract the results. (See ref. 2 for examples of this 
technique.) An alternative method for matrix prob- 
lems is described in appendix A of this paper. A 
method for estimating e A also is given in the ap- 
pendix of reference 5. 

Development of Optimum Step Size Formula 

A bound e on the total error is the sum of the 
truncation and condition errors 

-twi + s- (4) 

where \f!'\ is a bound on the second derivative in the 
interval [x,x+h\. In many cases, such a bound on 
the second derivative is not available, and a central 
finite difference approximation 

A a ^ = (5a) 

h s 

is used, so that 

e=\\*\ + \* A (5b) 

It should be noted that equation (5a) is an approx- 
imation to. f"(x) and is a reasonable bound on f" 
in the interval [x, x + h) only if f" varies slowly 
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near x. The step size h s used for estimating the sec- 
ond derivative is not necessarily the same as that 
used for estimating the first derivative. The optimum 
step size is found by minimizing the error defined in 
equation (4). Taking the derivative of e with respect 
to h, substituting $ for f", and setting the resulting 
expression equal to zero gives 


de 

dh 


M 

2 



Solving for h yields 


h 


opt 



( 6 ) 


(7) 


Clearly equation (7) is not applicable to deriva- 
tives of functions at points where the second deriva- 
tive is zero or close to zero — for example, when the 
function is nearly constant, linear, or an odd func- 
tion. For these special cases, the FD algorithm fails 
to provide a satisfactory finite difference step size. As 
discussed in reference 2, a computerized implementa- 
tion of the algorithm terminates with an appropriate 
error message when |$| is less than a certain thresh- 
old value. Equation (6) leads to the observation that 


h 


opt 
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( 8 ) 


Therefore, the optimum step size is the one which 
balances the truncation and condition errors. 


equation (10) for $ in equation (7). This results in 


ft. = 2(i + w)^HT (a) 

The interval h s is used to compute $ from equa- 
tion (5a); <J>, in turn, is used to calculate the er- 
ror c(<f>) from equation (9). If the error is too large 
(greater than 0.1), h s is increased. If it is too small 
(less than 0.001), h s is decreased. Once a value of h s 
is found which produces an acceptable value for $, 
the value of $ is used in equation (7) to produce an 
optimum step size, h op t- The limit of 0.1 is chosen to 
prevent the condition error from being excessive. The 
limit of 0.001 is chosen because when the condition 
error is very small, there is a risk of high truncation 
errors. 

Applications of FD Algorithm 

Sine Function 

A simple example to demonstrate the FD algo- 
rithm is sin x. This application is described in a 
tutorial fashion in appendix B to illustrate the steps 
in the procedure. 

Polynomial-Fit (Vandermonde) Problem 

The polynomial-fit example is based on the prob- 
lem of passing an iVth-order polynomial through 
N + 1 data points, Xj = i, t/ t - = /(xj), 

i = 0, 1, 2, ..., N. The form of the polynomial is 


Estimation of Second Derivative 

Before evaluating equation (7) for h op t, the sec- 
ond derivative is evaluated from equation (5a). This 
requires a value for h s , which is obtained by the fol- 
lowing trial-and-error procedure. The relative condi- 
tion error in f" resulting from use of equation (5a) 
is shown in references 1 and 2 to be bounded by 

am = 1*1 (9) 

The choice of an initial guess for h s is based on the 
conditions that the function is well scaled and that 
the function and its second derivative are comparable 
in magnitude, such that 


1 + 1/1 
(i + M ) 2 


( 10 ) 


where |x| is of the order unity. The trial value of h s to 
be used in equation (5a) is calculated by substituting 


N 

f( x ) = Y, ^ 

i=0 

To solve for the coefficients a,-, a system of N + 1 
equations in N + 1 unknowns is generated. It is 
required to calculate the derivatives of the coefficients 
dj with respect to the data point locations Xj. The 
equation for computing the coefficients is 


"1 

zo 

T 2 

x 0 

X 3 

x 0 ... 

t N i 

x 0 


' do ’ 


' 2/0 ' 

1 

Xl 

Ji 

X 1 

x 3 

x l 

r N 

X 1 


di 


2/1 

1 

x 2 

X 2 

x 2 

X 3 

X 2 . . . 

x N 

x 2 


d 2 


2/2 

1 

x 3 

T 2 

x 3 

X 3 

X 3 . . . 

x N 

X 3 

< 

d3 

• = . 

V3 

.1 

X N 

T 2 

X N 

T 3 

X N ... 

X N 

X N 1 


. a N . 


. 2 IN . 


( 12 ) 


The matrix in equation (12) is referred to as the 
Vandermonde matrix and is useful in the present 
work because the range of values (from unity to N n ) 
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causes it to become ill-conditioned even for relatively 
small values of N. For example, N = 10 generates 
values from 1 to 10 19 . 

In this example, the value for N was selected to be 
12 because it produced condition errors which were 
significant but not catastrophic. The values of xq 
and 1/0 were set to zero, so that ao = 0. The FD 
algorithm was used to find an optimum step size to 
use in approximating the derivatives of each coeffi- 
cient of the polynomial with respect to the location 
of the third data point, i.e., da^/dx 3, i = 1, 2, ..., 12. 
The data points ?/,• were generated from the equation 
2/j = i + i 5 + t 9 at values of i ranging from 0 through 
12. This method of generating data gives a means for 
estimating e a because the polynomial which exactly 
fits these data has coefficients a\ = 05 = ag = 1.0, 
and the remaining coefficients should be 0.0. Any 
deviation in the values of the coefficients is the result 
of numerical imprecision. A standard matrix subrou- 
tine denoted MATOPS was used to solve the system 
of equations. MATOPS offers three options (full piv- 
oting, partial pivoting, and no pivoting) correspond- 
ing to three levels of accuracy. The least accurate 
option was used for evaluating derivatives, and a first 
estimate for was obtained by subtracting the least 
accurate solutions from the most accurate. The sec- 
ond method used for estimating e a was the iterative 
improvement technique described in appendix A. 

The two estimates of e a were very close, but both 
were sensitive to the values of Xj. For the nominal 
point where Zj = i, e^ was evaluated as 8.3 X 10 -4 . 
However, at the nearby point of z,- = 0.01 -ft, (-A was 
estimated as 0.153. The low value of the error for 
Zj = i may be explained by a reduction in round-off 
errors when all the terms in the Vandermonde matrix 
are exact integers. 

Based on an value of 0.153, equation (11) gave 
h 3 = 2.2, and c(4) = 6.8 X 10 11 from equation (9). 
Because c(4>) was outside the interval [0.001,0.1], h s 
was reduced, successively, by factors of 10 till a value 
of 2.2 X 10 -3 resulted in c(4>) = 0.06. The optimum 
step size, h op t , was found from equation (7) to be 
1.71 x 10~ 4 . The derivative da\/dx 3 was computed 
by repeating the solution of equation (12) with Z3 
replaced by Z3 + h Q p t, while the remaining Zj values 
as well as the y t - values were left unchanged from 
their nominal values. The value of the derivative was 
obtained as —4.36028 X 10 6 . The exact value of the 
derivative was calculated analytically by Gaylen A. 
Thurston of NASA Langley Research Center to be 
—4.36003 X 10 6 , so that the error in the derivative 
at h op t is 250, or 0.0057 percent. The derivatives as 
a function of the step size are listed in table I for a 
wide range of h values bounding /i op t- It is evident 
that the choice of h op t is indeed a good one. 


TABLE I. FINITE DIFFERENCE EVALUATION OF 
d ai /dx 3 FOR POLYNOMIAL-FIT PROBLEM 


Step size, 
h 

Finite difference 
derivative 

Derivative 

•error 3 

Derivative 

error, 

percent 

10- 1 

-5.59003 x 10 6 

1.23 x 10 6 

28.2108 

10“ 2 

-4.46303 

1.03 X 10 5 

2.3624 

10 -3 

-4.37003 

1.00 X 10 4 

0.2294 

1.71 x 10- 4 (/.opt) 

-4.36028 

0.25 X 10 3 

0.0057 

10“ 4 

-4.35313 

0.69 x 10 4 

0.1583 

10“ 5 

-4.36095 

0.92 X 10 3 

0.0211 

10- 6 

-4.34963 

1.04 X 10 4 

0.2385 

10- 7 

-4.32533 

3.47 x 10 4 

0.7959 

10“ 8 

-4.05803 

3.02 X 10 5 

6.9266 

0 

1 

0 

i-H 

-5.31903 

9.59 X 10 5 

21.9953 


“Relative to exact value (—4.36003 x 10®). 

It was of interest to evaluate the error estimate 
of equation (5b) for this problem because the exact 
error was available. Using equations (5b) and (8) 
gives the bound for the minimum error e m ; n as 



So for = 0.153 and h op t = 1.71 x 10~ 4 ,e m j n = 
3580. This value is fourteen times larger than the ac- 
tual error. To resolve this large difference between 
the error predicted by the FD algorithm and the ac- 
tual error, the derivative was calculated in a nar- 
row range of step sizes in the neighborhood of h 0 p t. 
A sample of the results is given in table II. It was 
clear that near h op t, the value of the finite differ- 
ence derivative exhibits a random oscillation, and the 
small error obtained at h op t is somewhat fortuitous. 

The FD algorithm also was used for derivatives 
of the remaining polynomial coefficients with respect 
to Z3. Table III lists the second derivative step size 
h s , the optimum step size h 0 p t, and the value of 
for all 12 derivatives. It is seen from table III that 
hopt varies from 1.642 x 10~ 4 to 3.054 x 10~ 4 . 

Swept- Wing Problem 

A finite element structural model of a swept wing 
is shown in figure 1. The Engineering Analysis Lan- 
guage (EAL) finite element program (ref. 6) was used 
to model the wing. The model contains 194 elements 
consisting of rods, triangular membranes, and shear 
panels connected at 88 nodes and constrained at the 
root. Additional details of the wing model are avail- 
able in references 4 and 7. 
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TABLE II. FINITE DIFFERENCE EVALUATION 
OF dai/dx 3 IN THE NEIGHBORHOOD OF 
h op t = 1.71 x 10 “ 4 FOR POLYNOMIAL- 
FIT PROBLEM 


Step size, 
h 

Finite difference 
derivative 

Derivative 

error 0 

Derivative 

error, 

percent 

0.4356 X 10 -4 

-4.35986 x 10® 

170 

0.0039 

0.5445 

-4.35980 

230 

0.0053 

0.6806 

-4.35911 

920 

0.0211 

0.8507 

-4.36278 

2750 

0.0631 

1.063 

-4.35696 

3070 

0.0704 

1.329 

-4.35715 

2880 

0.0661 

1.662 

-4.36103 

1000 

0.0229 

1.714 (hopt) 

-4.36028 

250 

0.0057 

2.077 

-4.36072 

690 

0.0158 

2.596 

-4.36110 

1070 

0.0245 

3.245 

-4.36372 

3690 

0.0846 

4.056 

-4.36583 

5800 

0.1330 

5.071 

-4.36493 

4900 

0.1124 

6.338 

-4.36316 

3130 

0.0718 


“Relative to exact value (—4.36003 X 10®). 


TABLE III. VALUES OF t A ,h s , AND h opt FOR 
DERIVATIVES dajd 13 FOR POLYNOMIAL- 
FIT PROBLEM 


Coefficient 

£ A 

h s 

^opt 

“1 

0.1526 

2.209 x 10 -3 

1.714 x 10 “ 4 

“2 

0.4560 

5.395 

1.794 

“3 

0.5643 

0.600 

1.642 

“4 

0.3886 

4.982 

2.004 

“5 

0.1674 

2.314 

2.093 

a 6 

4.777 x 10 ~ 2 

1.748 

2.449 

“7 

9.273 x 10 -2 

0.770 

2.112 

“8 

1.231 x 10 -2 

2.807 

2.478 

a 9 

1.099 x 10 -3 

5.931 

2.884 

“10 

6.311 x 10 ~ 5 

2.010 

2.633 

“11 

2.102 x 10“ 6 

3.668 

2.913 

“12 

3.086 x 10-® 

4.444 

3.054 


The problem, as described in reference 4, con- 
cerns derivatives of the wing displacements with re- 
spect to element cross-sectional dimensions. In order 
to achieve accurate results, in reference 4 a time- 
consuming trial-and-error process was used to deter- 
mine an appropriate step size. The derivative of the 
normal displacement at node 44 with respect to one 
of the design variables is chosen to investigate the 
application of the FD algorithm. The design vari- 
able selected controls the thickness of the triangular 


membrane elements 1 through 6 located at the root 
of the wing model as shown in figure 1. 

The first step is to calculate an initial value for h s 
from equation (11). This calculation requires values 
of the function / and the independent variable x and 
an estimate of the error bound e A for this problem. 
In this example, x is the nominal value of the design 
variable (0.2 in.), / is the normal displacement of 
node 44 (40.8 in.), and e A (obtained by the method 
described in appendix A) is 1.0 X 10 -8 in. These 
values result in an initial h s value of 3.775 x 10 -5 in. 
used in equation (5a) to approximate <f>, which was 
then used to calculate c(4>) from equation (9). The 
value for c(<f>) was 0.06 and fell within the acceptable 
range [0.001,0.1]. The value of the optimum step 
size h 0 ptj 9-4 x 10 -6 in., produced a 0.0059 percent 
error in the finite difference derivative relative to 
the exact value obtained in reference 4. This step 
is better than the one (2 x 10 -5 in.) obtained by 
trial and error in reference 4 and was achieved with 
far fewer function evaluations. The variation of the 
error with step size is summarized in table IV and 
indicates that indeed the value of 9.4 x 10 -6 in. is 
close to optimal. Unlike the case of the polynomial- 
fit problem, for this problem equation (13) predicts a 
minimum error which is very close to the actual error 
(0.004 in. compared with 0.003 in.). 

TABLE IV. VARIATION OF DERIVATIVE OF TIP 
DISPLACEMENT WITH STEP SIZE FOR 
SWEPT-WING PROBLEM 


Step size,“ 
h, in. 

Error® in 
finite difference 
derivative, percent 

2 x 10 -3 

0.9400 

2 x 10“ 4 

0.0960 

c 2 X 10 -5 

0.0098 

d 9.4 x 10 _ ® 

0.0059 

2 x 10 ~® 

0.0200 


“Nominal thickness = 0.2 in. 

^Relative to analytically calculated value of —50.907 in. 
(ref. 4). 

“Best value obtained by trial and error in reference 4. 
d Step size calculated by FD algorithm. 


Extension of FD Algorithm to Vector 
Derivatives 

When the derivative of a vector with respect to 
a design parameter is required, the FD algorithm 
is not always useful in its original form. As seen 
in the polynomial- fit example (table III), the FD 
algorithm gives a different optimum step size for the 
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Figure 1. Geometry, element numbers, and node numbers for swept wing. Dimensions in inches. 
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derivative of each component of the vector. From 
the standpoint of efficiency, it is desirable in the case 
of a derivative of a vector to find a single step size 
for all components which is the optimum compromise 
among them. It is possible, of course, that no such 
compromise exists if the components have widely 
different condition and truncation errors. However, 
in most practical cases it is expected that such a 
compromise will be possible, and an extension of the 
FD algorithm described below is appropriate. 

In seeking a compromise among the step size re- 
quirements of the different components, it is natu- 
ral to seek to minimize some norm of the error vec- 
tor. Given an error vector E with components e t -, 
i = 1, 2, ..., N, there are three norms which may be 
used; 

IIEIIi = £ M (14) 

1 = 1 


to the polynomial-fit problem discussed earlier pro- 
duced 12 optimum step sizes, one for each coefficient. 
(See table III.) To calculate all 12 finite difference 
sensitivity derivatives with the same step size, the 
extended FD algorithm was used. To apply the ex- 
tended FD algorithm, the second derivative step size 
h s obtained for ai (2.209 X 10~ 3 , see table III) was 
used for all 12 coefficients. All 12 c(<&) values were 
in the acceptable range of [0.001,0.1]. The values 
for h 0 pt changed slightly from the values given in ta- 
ble III and ranged from 1.714 x 10 -4 to 3.364 x 10 -4 . 
Subdividing this range into 20 intervals resulted in 
21 vectors, each containing 12 error components. The 
IjEHj and ||E|| 2 norms (eqs. (14) and (15)) produced 
the same step size of 1.879 X 10 -4 , whereas the ||E[| ^ 
produced a step size of 1.797 X 10 -4 . A check of the 
actual errors based on these step sizes showed that 
all 12 errors are acceptably small. 


||E|| 2 = VNe ims = 


N 



(15) 


||E||oo = max |e,-| (16) 

t 

where e rms denotes the root-mean-square of the 
error. 

The extended FD algorithm proceeds through the 
following steps: 

1. Select an initial h s (e.g., by using eq. (11)), and 
calculate $ from equation (5a) for each component 
of the vector. 

2. Estimate the condition error in $ (c(<I>) from 
eq. (9)) for each component, and vary h s to bring 
that error into the required interval of [0.001,0.1]. 
Failure to find a common h s is an indication that it 
is futile to use a single step size for all components. 
In that case, the components may be divided into 
two or more groups and subsequent steps repeated 
for each group separately. 

3. Estimate the optimum step size /i 0 pt for each 
component from equation (7). The largest and 
smallest values of h op t are denoted h max and h m ; n , 
respectively. 

4. Divide the interval [h min , h max ] into m equal 
subintervals resulting in m + 1 step size values h\ = 
^mim h 2, ^3, •••> ^m+l = V ax- Equation (4) is used 
to estimate the error in each component for each /i t - 
value and the associated values of $ and e A . 

5. Select as h op t that value of h{ which minimizes 
the desired error norm. 


Applications of Extended FD Algorithm 


Polynomial-Fit Problem 

The application of the standard FD algorithm 


Swept-Wing Problem 

To demonstrate the application of the extended 
FD algorithm to the swept wing, a set of 12 displace- 
ment derivatives were chosen. These are the deriva- 
tives of the normal displacements at nodes 5 through 
10 and 39 through 44. (See fig. 1.) The FD algo- 
rithm was first used to find 12 optimum step sizes. 
For this purpose, h s = 3.775 x 10~ 5 in. was used for 
all 12 components, and one step size was extracted 
for calculating all 12 derivatives. Unlike the previous 
case, the three norm criteria produced different step 
sizes: 9.05 x 10 -6 in. from ||E||i, 9.44 x 10 -6 in. from 
||E|| 2 , and 9.6 x 10 -6 in. from ||E||oo. Derivatives 
using these step sizes were compared with analytical 
results from reference 4. All three gave sufficiently 
accurate results. For example, table V shows that a 
step size of 9.44 X 10~ 6 in. produced excellent results 
compared with the analytical solution. 

Effect of e A 

Because the error bound c A may be uncertain in 
some problems, the effect of e A on the step size cal- 
culations and derivative results is of interest. To test 
this effect, the value of the nominal e A was varied by 
a factor of 10 for each of the single-derivative prob- 
lems. The results are shown in tables VI, VII, and 
VIII for the sine function, the polynomial-fit prob- 
lem, and the swept-wing problem, respectively. In 
each case, the modified values of e A had a noticeable 
but small effect on the results, and the errors were 
still acceptable. 

Concluding Remarks 

This paper deals with methods for obtaining near- 
optimum step sizes for finite difference approxima- 
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TABLE V. FINITE DIFFERENCE DERIVATIVES 
OF SWEPT-WING TIP DEFLECTION FOR 
/i op t = 9.44 x 10 — 6 in. 


Node 

number 

Derivative 

Error, 

percent 


Analytical method 
of reference 4 

5 

-2.9372 

-2.9374 

0.00681 

6 

-3.7578 

-3.7580 

0.00532 

7 

-4.5608 

-4.561 

0.00438 

8 

-5.5679 

-5.5682 

0.00538 

9 

-10.645 

-10.646 

0.00939 

10 

-11.597 

-11.598 

0.00862 

39 

-47.737 

-47.739 

0.00419 

40 

-48.219 

-48.222 

0.00622 

41 

-49.773 

-49.776 

0.00603 

42 

-50.112 

-50.115 

0.00599 

43 

-50.496 

-50.499 

0.00594 

44 

-50.904 

-50.907 

0.00589 


TABLE VI. EFFECT OF e A ON DERIVATIVE OF sin x 



^opt 

Finite 

difference 

derivative 

Derivative 
error, a 
percent 

Nominal 
(2.188 x 10 -7 ) 

1.108 x 10 -3 

0.706679 

0.06 

0.1 x Nominal 
(2.188 x 10“ 8 ) 

1.20 x 10 -4 

0.708333 

0.17 

10 x Nominal 
(2.188 X 10 -6 ) 

3.464 x 10“ 3 

0.705831 

0.18 


“Relative to exact value (0.707107). 


TABLE VII. EFFECT OF e A ON DERIVATIVE 
dax/dx 3 IN POLYNOMIAL-FIT PROBLEM 


C A 

^opt 

Finite 

difference 

derivative 

Derivative 

error, 0 

percent 

Nominal 

(0.153) 

1.714 X 10~ 4 

-4.3603 X 10® 

0.006 

0.1 X Nominal 
(0.0153) 

2.656 x 10-® 

-4.3596 X 10® 

0.009 

10 x Nominal 
(1.53) 

4.885 x 10~ 5 

4.3636 x 10® 

0.08 


“Relative to exact value (—4.3600 x 10®). 


tions to first derivatives with particular application 
to sensitivity analysis. A technique denoted the fi- 


TABLE VIII. EFFECT OF e A ON DERIVATIVE IN 
SWEPT-WING PROBLEM 


«A> in - 

h Q pt> 

Finite 

difference 

derivative 

Derivative 

error, 0 

percent 

Nominal 

(0.10348 x 10 -7 ) 

9.4 x 10 -6 

-50.904 

0.0059 

0.1 X Nominal 
(0.10348 X 10~ 8 ) 

2.525 x 10“® 

-50.914 

0.0138 

10 X Nominal 
(0.10348 x 10~®) 

2.929 X 10 -5 

-50.900 

0.0138 


“Relative to exact value (—50.907). 

nite difference (FD) algorithm, previously described 
in the literature, is reviewed and applied in a tuto- 
rial manner to the derivative of a sine function. The 
original FD algorithm is applicable to one derivative 
at a time. This paper describes an extension of the 
FD algorithm to the calculation of several derivatives 
of matrix equations. 

The FD algorithm requires an estimate of the er- 
ror in the computed value of the functions being dif- 
ferentiated. This paper describes a simple method 
for calculating the error estimate for matrix prob- 
lems. Denoted the iterative improvement technique, 
the method involves a second solution of the matrix 
equation. Since it uses the same decomposed matrix 
that was used in the original solution, it imposes a 
minimal computational burden on the analysis. 

Both the original and extended FD algorithms 
were applied to sensitivity analysis for a data-fitting 
problem in which a polynomial is passed through 
data points, and the derivatives of the coefficients 
of the polynomial with respect to uncertainties in 
the data were calculated. The algorithms also were 
applied to sensitivity analysis of the structural re- 
sponse of a finite-element-modeled swept wing. In a 
previous paper, this sensitivity analysis of the swept 
wing required a time-consuming trial-and-error effort 
to obtain a suitable step size, but it proved to be a 
routine application for the FD algorithm herein. 

One limitation of the FD algorithm is that it re- 
quires a finite lower bound on the second derivative 
of the function being differentiated. This bound is 
generally estimated by a central difference approxi- 
mation. In certain special cases, such as when the 
function is nearly constant, linear, or an odd func- 
tion, the second derivative is zero or close to zero, 
and in such a case, the FD algorithm is not a reli- 
able provider of an optimum step size for the first 
derivative. In the absence of this special condition, 
however, the FD algorithm has been found to be use- 
ful, accurate, and reliable. 


8 










Appendix A 

Iterative Improvement Technique for 
Estimating e A 

A method for estimating e A is based on the so- 
called “iterative improvement” of a solution to a 
system of linear equations. This process is explained 
for the system 

KU = F (Al) 

where K is an n x n matrix, and F and U are 
n x 1 vectors. In finite element structural analysis, 
K is the stiffness matrix, U is a displacement vector, 


and F is a force vector. If U is the solution obtained 
numerically, F can be defined as 

F = KU (A2) 

Subtracting equation (A2) from equation (Al) gives 

K(U - U) = F - F (A3) 

so that an estimate of U — U may be obtained by 
solving equation (A3). The solution gives an order 
of magnitude estimate for the error e^, and since it 
uses the same decomposed matrix that was used in 
the original solution (eq. (Al)), it imposes a minimal 
computational burden on the analysis. 
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Appendix B 

Application of FD Algorithm to Sine 
Function 

A simple example to demonstrate the FD algo- 
rithm is sin x. The derivative of this function was 
evaluated with a hand calculator which carries up to 
10 significant digits. The function sin x is convenient 
because the exact derivatives are known and can be 
compared with the approximations. The bound on 
the absolute error in the computed function values e A 
is controlled by the number of digits used on the cal- 
culator. In this example, the number of digits used 
was six. The function sin x and its first derivative 
were evaluated at x = 0.785398 rad, or 45°. There- 
fore, /( x) = sin x = 0.7071067812 to 10 significant 
figures, and sin x = 0.707107 to 6 significant fig- 
ures. The difference between these two numbers is 
e A = 2.188 x 10 -7 . 

The first step is the calculation of the initial step 
size h s to be used in estimating the second derivative. 
From equation (11) 


' , ’ = 2 <i + |l l>V / T?7f 


Since c(3>) is outside the acceptable interval 
[0.001,0.1], h s must be increased, and new values for 
$ and c($) must be calculated. For convenience, h s 
was increased by a factor of 10. The new value for 
|3>| is 0.711656 for c($) = 0.008621. Since this value 
falls within the acceptable interval, may be used 
to calculate h Q pt from equation (7). 



= 2 


2.188 x 10-7 
0.711656 


= 1.108 x 10~ 3 


Substituting h op t into the forward difference for- 
mula (eq. (1)) gives an approximation to the first 
derivative 


f\ x ) ~ /(* + ft opt) ~ /(») 
*opt 

_ 0.707890 - 0.707107 
1.108 x 10 -3 
= 0.706675 


= 2(1 + 0.785398) 


/ 2.188 x lO- 7 
1 + 0.707107 


= 0.001278 


Using h s in equation (5a) gives 


Comparing this result with the exact derivative in- 
dicates that the use of h 0 p t gives a relative error of 
0.06 percent. Table BI shows derivatives calculated 
with step sizes larger and smaller than the optimum 
along with the error relative to the exact derivative. 


_ f(x + h s ) - 2 /( x) + f(x - h s ) 

0.708010 - 2(0.707107) + 0.706202 
“ 0.001278 2 

= -1.225 


The relative condition error is calculated from equa- 
tion (9) 


*(*) = 


4(2.188 X lO -7 ) 
(0.001278) 2 (1.225) 

0.4374 


TABLE BI. FINITE DIFFERENCE DERIVATIVES OF 
sin x AT 45° FOR VARIOUS STEP SIZES 


Step size, 
h 

Finite difference 
derivative 

Derivative 
error, 0 percent 

1.108 X 10 _1 

0.666525 

5.74 

1.108 X 10 -2 

0.703159 

0.56 

1.108 x 10“ 3 (/.opt) 

0.706675 

0.06 

1.108 X 10~ 4 

0.703971 

0.40 

1.108 X 10~ 5 

0.722022 

2.11 


“Relative to exact value (0.707107). 
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